clc; clear; close all;

alpha = 2.4;
A = [-alpha, 0, 0;
      0,     0, 0;
      0,     0, 0];
B = [1, -4, -3.5;
     0,  1,  2;
    -1, -4,  1.5];

k_list = linspace(0, 0.2, 400);
dt = 0.01;
t_total = 600;
transient = 550;
trans_steps = round(transient/dt);

init_state = [1e-6; 0; 0; 0; 0; 0];
init_inner = [0; 0; 0];  % q11, q21, phi0=0

% 预分配cell数组存峰值和对应k
phi_peaks_cell = cell(length(k_list), 1);
k_phi_cell = cell(length(k_list), 1);

q11_peaks_cell = cell(length(k_list), 1);
k_q11_cell = cell(length(k_list), 1);

options = odeset('RelTol',1e-8,'AbsTol',1e-10);

parfor idx = 1:length(k_list)
    k = k_list(idx);
    x0 = [init_state; init_inner];
    tspan = 0:dt:t_total;
    [t, Y] = ode45(@(t,x) mCNN_coupled_ode(t,x,A,B,k), tspan, x0, options);
    
    % 去暂态数据索引
    idx_ss = t > transient;
    
    phi = Y(idx_ss, 9);
    q11 = Y(idx_ss, 7);
    
    pks_phi = findpeaks(phi);
    pks_q11 = findpeaks(q11);
    
    if ~isempty(pks_phi)
        phi_peaks_cell{idx} = pks_phi(:);
        k_phi_cell{idx} = repmat(k, length(pks_phi), 1);
    else
        phi_peaks_cell{idx} = [];
        k_phi_cell{idx} = [];
    end
    
    if ~isempty(pks_q11)
        q11_peaks_cell{idx} = pks_q11(:);
        k_q11_cell{idx} = repmat(k, length(pks_q11), 1);
    else
        q11_peaks_cell{idx} = [];
        k_q11_cell{idx} = [];
    end
end

% 合并所有点
phi_peaks_all = vertcat(phi_peaks_cell{:});
k_phi_all = vertcat(k_phi_cell{:});

q11_peaks_all = vertcat(q11_peaks_cell{:});
k_q11_all = vertcat(k_q11_cell{:});

% 绘图
figure; hold on;
scatter(k_phi_all, phi_peaks_all, 4, 'b', 'filled');
scatter(k_q11_all, q11_peaks_all, 4, 'r', 'filled');
xlabel('k');
ylabel('q_{11} & \phi');
legend({'\phi', 'q_{11}'});
grid on; box on; hold off;

% 系统微分方程
function dxdt = mCNN_coupled_ode(~, x, A, B, k)
    X1 = x(1:3);
    X2 = x(4:6);
    q11 = x(7);
    q21 = x(8);
    phi = x(9);

    M1 = exp(-q11);
    M2 = exp(-q21);
    Y1 = M1 .* X1;
    Y2 = M2 .* X2;

    W_phi = 1 - 2*exp(phi);
    K = diag([0, 0, k]);

    dX1 = -X1 + A*Y1 + B*X1 + K*W_phi*(X1 - X2);
    dX2 = -X2 + A*Y2 + B*X2 - K*W_phi*(X1 - X2);
    dq11 = X1(1);
    dq21 = X2(1);
    dphi = X1(3) - X2(3);

    dxdt = [dX1; dX2; dq11; dq21; dphi];
end

